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ABSTRACT 


This study adapts some established target tracking techniques for use in the 
maritime surface surveillance role and tests them with computer generated data. Computer 
simulation of a track splitting tracker capable of operating in this undersampled and 
asynchronous environment is presented. The tracker uses standard and extended Kalman 
Filter algorithms to estimate target state from latitude and longitude or line of bearing 
position measurements. Prior to state estimation, all measurements are processed to retain 
only those that meet feature and geographic gate thresholds. All measurements passing 
these criteria will update the target state and be scored based on a goodness of fit with the 
model. The state estimate with the best score is selected as the correct one for display 
purposes, while all state estimates continue to be processed with subsequent measurements. 
Several runs of the simulation are discussed here to illustrate the performance of track 


splitting and the effect of several key tracker parameters. 
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I. INTRODUCTION 


A. BACKGROUND: THE MARITIME SURVEILLANCE ENVIRONMENT 

The complexity of the target tracking problem varies according to a number of 
factors. Significant research has been devoted to the study of trackers in varying levels of 
target motion model, (plant) noise and measurement noise. Performance under a variety of 
target maneuvers has also been studied. Clutter environments have been added to the 
model, further complicating the tracking problem. Tracking problems containing multiple 
targets, multiple sensors, and target maneuvers are most complex, and extensive research 
has gone into developing algorithms for these situations. A maritime surveillance 
application for such a tracker is the focus of this study. The characteristics of this 
environment are unique in many respects, presenting an interesting and challenging tracking 
problem. 

One of the most problematic aspects of maritime surveillance tracking is the long 
time interval between position updates. Typical air target tracking systems, such as a fire 
control system, will receive measurements of the target many times each second from 
tracking radars with high pulse repetition rates. A maritime surveillance tracking system, 
however, will receive measurements of contacts at intervals measured in minutes or hours 
depending on such factors as contact geographic location, asset availability, and source level 
or signal strength. Fortunately, because it is a surveillance system, the maritime 


surveillance tracker has a relaxed accuracy requirement when compared to an air target 


tracker: accuracy within hundreds or a few thousand yards is acceptable in this application, 
but clearly unacceptable in air target tracking. Also, the relatively slow speed of a surface 
contact mitigates the effect of longer update intervals. Nonetheless, the undersampled 


environment is a significant challenge that a maritime surveillance tracker must overcome. 


There are a variety of participants in the maritime surveillance effort: active radars 
from ships or aircraft at sea, nationally tasked sensors, shipping control officials using 
civilian shipping reports, to list some. While the sources or sensors from which a contact’s 
position is obtained may be dissimilar (acoustic or electronic line of bearing, active radar, or 
a filed plan of intended movement), the data is integrated into the maritime surveillance 
picture via the Sensor Report. From this report, the tracker receives either a line of bearing 
or latitude and longitude for the target and a confidence region. This confidence region may 
be an ellipse (for latitude and longitude reports) or a bearing swath (for line of bearing 
reports) and is sensor dependent. Other information that may be contained in the sensor 
report includes estimated course and speed and feature data such as emitter frequency and 
pulse repetition frequency. Sensor reports allow dissimilar sensors to participate in the 


tracking effort. 


Unlike many tracking problems where position information is received at regular 
intervals, there is no synchronization of sensor reports in maritime surveillance. The 
various sensors participating each operate according to their own reporting abilities and 
frequently have irregular reporting periods. | Consequently, batch processing of 


measurements is impractical. 


Finally, the number of tracks that a maritime surveillance tracker must handle to be 
of value further complicates the design. With the goal of providing situational awareness of 
the surface traffic in a geographic area to naval forces, a maritime surveillance tracker 
would have to support thousands of tracks. While advances in computing dnd memory 
storage technologies will enable the design of a high volume tracker, the number of contacts 
to be tracked in a geographic area will vary with time such that the tracking system is faced 
with an unknown number of targets. This restriction eliminates many tracking algorithms 


from consideration. 


B. THE TRACK SPLIT APPROACH 

The track split approach is a multi-target tracking algorithm that lends itself to 
several problems, including single target in clutter, multiple targets in clutter, and multiple 
targets in clutter with the number of targets unknown. Unlike some simple trackers which 
assume the measurement closest to the estimate is the correct one, track splitting trackers 
will use all measurements that fall within a validation region, or gate, to update the target 
state. A likelihood function is then used to eliminate unlikely sequences that presumably 
represent false tracks. In this way, the tracker postpones important decisions about which 
measurements originated from a target until additional measurement information is 
received. 

The track split approach shares the Kalman Filter based prediction and update 
sequence used by other trackers. This widely used technique obtains the optimal estimate of 


the target state, for a given dynamic model, by recursively minimizing the mean square 


error. ‘The filter is tunable to a variety of applications and environments by properly 
selecting the dynamic and noise model. Asynchronous or missed measurements are readily 
handled by the filter and a measure of tracker accuracy (relative to the specified model) is 
always available in the form of the target’s covariance matrix. Because of its versatility, a 
Kalman Filter based Track Splitting tracker can be developed for a maritime surveillance 
role. 

The track split approach was selected for this study because it is the simplest 
algorithm capable of meeting the requirements imposed by the maritime surveillance 
environment. Track splitting is a technique for allowing the Kalman Filter to handle false 
measurements and multiple contacts. Furthermore, this approach does not require the 
number of targets to be known and will automatically initialize new tracks when more than 
one measurement falls within the gate of a currently held track. Although other algorithms 
are also capable of performing under these conditions (most notably, the Multiple 
Hypothesis Tracker (MHT)), the track split approach is a logical starting point in the design 


and evaluation of a maritime surveillance tracker. 


e: THESIS OUTLINE 

The remainder of this thesis is organized as follows: Chapter If will discuss the 
contact and measurement model by presenting the state space dynamic representation and 
noise model. Chapter II will detail the tracker algorithm, including a review of standard 
and extended Kalman Filter equations, track split logic, gating, and estimate error ellipses. 


Chapter IV will present and discuss the simulation parameters and results. Chapter V will 


draw conclusions about this study and outline some areas for future research. The tracker 


simulation code, written for MATLAB®,, is contained in the Appendix. 





II. CONTACT AND MEASUREMENT MODEL 


A. CONTACT MOTION MODEL 


The tracker assumes contacts have a two-dimensional (x-y), second order (position 
and velocity) dynamic model. The model is implemented in discrete form according to: 
ee wee Re, Coe ae Vay Pek 
where x is the state space vector composed of x position, x velocity, y position, and y 


velocity. F is the linear matrix: 


ee a 
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Here, AT is the time interval between samples. Plant noise is included in the 
dynamic model to account for contact deviations from straight-line motion. Ux 1S a zero 


mean Gaussian random variable with known covariance, Q , given by: 
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The parameter g 18 used to adjust the magnitude of plant noise. Properly selecting 
this parameter allows the designer to optimize tracker performance for the expected 
behavior of the contacts. Contacts that maintain relatively straight line motion, in open 
ocean transit for example, will not deviate greatly from the dynamic model. A low value 
for g will give optimum tracker performance. Maneuvering contacts, however, conform to 


the model poorly. A large value for g° allows the tracker to maintain track. 


B. MEASUREMENT MODEL 

Sensors feeding into the maritime surveillance tracking system obtain contact 
positions as either single dimensional lines of bearing or two dimensional range and bearing 
measurements. Passive systems, such as acoustic hydrophones or electronic direction 
finding receivers, will produce lines of bearing. Active radars and electronic geolocation 
systems employing techniques such as Time Difference of Arrival measurements produce 
two-dimensional position reports. All sensors have an associated measurement error which 


is modeled as a zero mean Gaussian random variable. 


Position reports contain the sensor’s measurement, along with a metric of 
measurement uncertainty. In the case of Line of Bearing reports, a 95% confidence swath 
represents about 2.5 standard deviations of the measurement noise. Thus, 69 can be readily 
found from the sensor report. Latitude and Longitude reports will report measurement error 
statistics in the form of an error ellipse. This ellipse 1s related to the eigenvalues and 


eigenvectors of the sensor’s covariance matrix, R, (Kirk, 1975): 





+ ——=C wal 


where 
A, is the eigenvalue associated with the X’ eigenvector of R 
Az is the eigenvalue associated with the Y’ eigenvector of R 
c’ is aconstant that determines the size of the ellipse 


The relationship between c” and the probability of a measurement lying within the 
ellipse is tabulated for the two dimensional case (Kirk, 1975). In this study, c” is chosen to 
be 6, representing approximately 95% probability that the measurement falls within the 


error ellipse. 


The tracker must invert the ellipse generating process to obtain the measurement 
covariance matrix. The sensor report contains the semi-major axis length (a), semi-minor 


axis length (b), and orientation of the ellipse (¢). The eigenvalues can be found by: 


CS = 
ZG 
25 
b 
A, =— 
: ZG 


The orientation of the semi-major and semi-minor axes may be expressed as unit 


vectors, giving the eigenvectors of the matrix that generated the ellipse directly. A similar 


matrix to R (one with the same characteristic equation) can be obtained as follows (Kaplan, 


1981): 


Re BG. 2.6 
where C is composed of the eigenvectors obtained from the ellipse orientation: 
Cc=[x’ Y’ 2a 


and B is composed of the eigenvectors obtained from the semi-major and semi-minor axes: 


B=| uy 2.8 
0 4, 


While the matnx R’ is not necessarily identical to the sensor’s covariance matrix, it has the 
same characteristic equation and is therefore acceptable for use in the tracker. 

A coordinate transformation from polar (range and bearing) to cartesian (latitude 
and longitude) occurs at the sensor. For purposes of this simulation, the cartesian position 
in distance units is transformed into latitude and longitude coordinates using a flat earth 
model, with one degree equal to 60 nautical miles (nm). A nonlinear rotation is used to 


obtain the measurement in cartesian: 


xX ‘s r sin Osx 2.9 
Y TacOSs ay = 


where r,, and 6@,, are the measured range and bearing to the contact. 6,, is the angle measured 
from a vertical axes, or "True North" reference. The measurement covariance matrix must 
also be transformed from polar to cartesian coordinates. The resulting elements of the 


cartesian covariance matrix (Bar Shalom and Li, 1995): 
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7) ; 
R,, =(o0°-—r’o;)sin 6, cos 6, 


where Oo, and og are the range and bearing standard deviations of the sensor. For the 
simulation, the measurement and covariance are transformed from polar to cartesian 
coordinates by the m-file pole2cart. The listing is contained in the Appendix. 

Equipped with these models, a tracker can be implemented that will produce 
estimates from line of bearing or cartesian position reports. The measurement and sensor 
covariance can be determined either directly, as in the case of the line of bearing report, or 
by analyzing the error ellipse associated with a cartesian report. Since measurements are 
generated for the simulation in this thesis, a technique for linearizing range and bearing 
measurements to obtain the position in cartesian coordinates is also required. The 
development of a track splitting tracker in the upcoming chapter will draw on the models 


discussed here. 


Il 





It. THEORETICAL BASIS FOR THE TRACK SPLITTING TRACKER 


A. THE KALMAN FILTER 

The Kalman Filter is a recursive algorithm for estimating the contact’s current state 
vector from its past estimate and current noisy measurement. With accurately modeled 
plant dynamics and sensor noise, the algorithm minimizes the mean square error in the 
resulting estimates. Several excellent sources exist for detailed derivations of the Kalman 
Filter (Bar Shalom and Li, 1998). Only the resulting equations are presented here. The 


following definitions and notations will be used: 


X41, = the estimate of state vector (x) at time k + 1 given the measurement 


at time k (prediction) 
X,, = the estimate of state vector (x) at time k given the measurement 


at time k (corrected estimate) 
P 


k+ilk 


P.,. = the covariance of state vector (x) at time & given the measurement at time k 


= the covariance of state vector (x) at time k + 1 given the measurement at time k 


The Kalman Filter recursion contains a prediction and a correction step. The 
mechanics of these steps are now described. 

Prediction: The weighted least squares estimate and associated covariance matrix 
for the next sample time are predicted based on plant and noise models according to these 
equations: 


Xesye = FLX y 
ull 


P ein = FP ay F, +Q, 
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Correction: The estimate and covariance matrix at time k+/ are corrected with the 
new (time k+/) measurement. 


Xp eer = ik a K 412 k+l 
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The standard Kalman Filter presented above can be used when the plant and 
measurement models are linear, as in the case of cartesian position reports. The 


measurement matrix, H, for the standard filter is: 
10 0 0 
= 35 
00 1 0 


The line of bearing report, however, requires a nonlinear measurement matrix since 


O@(x) = arctan ts 3.4 


x; 
The Extended Kalman Filter provides a means to handle nonlinear state or measurement 
equations by linearizing the nonlinear equations about the estimate. For the line of bearing 
problem, the linearized measurement matrix obtained from a first order Taylor Series 


expansion Is: 


x — x 
Ee per as 0 3.5 
XxX, tans Xt x eee 
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The measurement matrix is linearized about the predicted state and must be calculated with 
each new measurement. Once the linearized measurement matrix is obtained, it may be 
used in the correction step in the same manner that the static measurement matrix (Eqn. 3.3) 
is used with the standard Kalman Filter to obtain the corrected estimate. 

Since the Kalman Filter is recursive, an initial state estimate and state covariance are 
required to begin the prediction and correction process. A first order polynomial curve fit 
method for obtaining these from two (cartesian) measurements may be used (Bar Shalom 
and Li, 1998). The latest measurement is used for the position estimate, while the velocity 


estimate is obtained from the difference in position. The initial covariance is calculated: 


l — 0 0 
] 0 0 0 AT 
= l 
P, = 
0 =) 
; am 2 As cart 1 0 nis 0 0 
AT ST ‘0 0 0 —| 


Rear is as derived in the previous chapter (Eqn. 2.10), with the subscript denoting the 
measurement with which the covariance matrix is associated. The m-file kal2init generates 
the filter initializations from the plant model and simulation parameters. 

The Kalman Filter recursions form the backbone of the tracker algorithm, predicting 
and correcting the contact’s state estimate based on the most recent state estimate, the plant 
model, and the most recent measurement of the contact. In the multitarget tracking 
problem, some data association decisions must be made to pair measurements with existing 
tracks before this process can begin. A technique for comparing the features of the 


measurement with the known features of the target can be employed to quickly reduce the 


IS 


number of candidates for update. Also, a gating procedure will eliminate the most unlikely 


candidates by considering only the measurements within a region around the prediction. 


B. MEASUREMENT REJECTION USING FEATURES 
One technique employed in this study relies on the assumption that the features for 
both a contact being tracked and a false measurement are normally distributed with different 


statistics. The vector of features used here is: 


Bh 


emitter frequency (GHz) 
emitter pulse repetition frequency (pps) 


Clearly, not all sensors would include these features in the position report, so clutter 
rejection using features can not always be accomplished. The statistics of the feature 


random variables are denoted: 


ae 


E|prf | 3.8 


Distances are then formed to compare each measurement with statistics of the contact and 
clutter models. Here, the subscnpt "7" denotes target, and "C" denotes clutter, or false 


target: 


Q 
ae | 
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Small distances indicate a measurement that closely matches the feature model, based on 
the feature statistics. Forming a ratio and establishing a threshold provides a criterion for 


rejecting a measurement as clutter: 


La > Threshold 3.10 





C 


Alternatively, when the clutter feature means are not known, a gating technique can 
be used with the target feature distance alone. In this case, measurements are rejected when 
dr exceeds a Chi-Squared distribution threshold. Using a threshold of 6 will retain the 
correct measurement with 95% confidence while discarding measurements with features 
that conform poorly with the expected values. With both feature rejection techniques, 
measurements that meet the feature criteria remain valid candidates for update. Further 
processing will reduce the number of candidates based on the state estimate and position 


measurement. 


c. ELLIPSOIDAL GATING 

Just as a confidence ellipse exists around a noisy measurement, a region can be 
defined around an estimate that will establish the area within which a candidate observation 
must lie with some probability. The smallest region for a given probability is an ellipse 
defined by the innovation and innovation covariance resulting from a measurement, (Bar 


Shalom and Li, 1995): 
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where Sah 
< = HP,,.,H’+ R = the innovation covariance 
y = chi - squared threshold 


Z—Z,,,, = the measuremen t’s innovation 


The above quadratic (norm of the innovation squared) is a Chi-Squared random variable 
with the number of degrees of freedom equal to the dimension of the measurement (two for 
cartesian, one for line of bearing). Thus, for each measurement received and each track 
held, the measurement innovation and innovation covariance must be calculated. With 
these, the norm of the innovation squared may be calculated and compared to the threshold. 
Norms exceeding the threshold indicate a pairing that is highly unlikely and need not be 
considered. 

The application of this gating quadratic is straightforward for the cartesian position 
report. In these simulations, a threshold value of 6 has been used for cartesian reports, 
representing a 95% confidence gate. For line of bearing reports, the linearized measurement 
matrix must be used in determining the innovation covariance. The 95% threshold for this 
single degree of freedom distribution has been approximated as 4 in this simulation, but in 


reality is slightly less. 


D. TRACK SPLITTING AND THE LIKELIHOOD FUNCTION 
If only one measurement remains a candidate for track updating after processing 


with feature measurements and gates, no further tracker logic is required. The candidate 
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measurement is passed off to the Kalman Filter (standard or extended) and the corrected 
estimate can be obtained. If no measurements remain a candidate after processing, options 
can be included in the tracker to increase thresholds in the feature or ellipsoidal gate 
algorithms in an effort to obtain a candidate. (No such logic has been included in this 
simulation, however). When more than one measurement remains a candidate, the tracker 
is forced to make a decision as to which one most likely originated from the contact being 
tracked. With the track split approach, this decision is deferred until additional information 
(measurements) can be obtained. The predicted state is updated with each measurement (by 
the Kalman Filter correction step) so that a single track will become n tracks after 
processing, where n is the number of candidate measurements. For each of these tracks, a 
Kalman Filter prediction will be obtained so that a gate may be established for subsequent 
measurements. Feature processing, gating, splitting, and updating will be repeated for each 


track as new measurements are received. 


With this logic, the number of tracks held in the system can grow exponentially. To 
aid in identifying the correct track, (that is, the one most closely matching the plant model) a 
likelihood function based on the Gaussian assumptions of the plant and measurement model 
is computed recursively each time a measurement updates a track (Bar Shalom and Li, 


1998): 
A, = 4.4 t lz i A pee lz nv zap 3.42 


That is, the likelihood that a sequence of measurements originated from a track conforming 


to the model is the sum of the previous score and the norm of the innovation squared. A 
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low score indicates that the sequence conforms to the model well and is likely a valid track. 


Thus, a metric exists to aid in managing the number of tracks in the system 


E. TRACK DROPPING CRITERIA 


A Track Dropping algorithm is essential for a track splitting tracker to prevent 
overloading of the computational and data storage capacities of the system. Track scores, 
computed by the likelihood function (Eqn. 3.12), are one available tool for selecting the 
state estimates that have arisen out of the most unlikely measurement sequences. By 
choosing an appropriate threshold, the number of false tracks held in the system can be kept 
to a minimum without inadvertantly dropping valid tracks. Another parameter that aids in 
track management is the most recent time at which an estimate was updated with a 
measurement, or update time. A threshold can be set for this parameter as well, triggering 
the tracker to automatically eliminate a held track once it is exceeded. A combination of 


these techniques is included in this simulation. 
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IV. COMPUTER SIMULATION 


A. OVERVIEW 

The listing for a simulation of a track splitting tracker written for MATLAB® is 
contained in the Appendix. This tracker incorporates the theoretical models previously 
discussed. Figure 1 shows the functional flow of the tracker simulation. Two dimensional 
constant velocity motion is simulated for two targets in a flat earth coordinate system. An 
instantaneous turn may be specified by selecting a time parameter for that occurrence. Two 
attributes, emitter frequency and pulse repetition frequency (PRF), may also be associated 
with the contacts. A range and bearing measurement of contact position from a fixed sensor 
position is made nearly every 90 minutes. The contact’s features are also obtained at each 
measurement. Bearing-only measurement times may also be specified. Along with the 
noisy measurements of contact position, false measurements that did not originate from the 
contact are generated at sample times in an area around the contact’s position. These false 
measurements have their own associated feature measurements. Thus, a position 
measurement and the associated feature measurements at a specific measurement time can 


be considered a single sensor report. 


After generating all measurements at a sample time, track processing begins. The 
feature rejection loop compares the expected feature values for a given track with the 
measured values for each report. Reports that do not conform well to the track’s expected 


features are removed from the set of candidates for the track being processed, but will be 
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Figure 1. Functional Flow of Track-Splitting Simulation 


considered for subsequent tracks. The remaining candidate reports are filtered through a 
geometric gate around the track’s expected position. For all measurements that fall within 
the gate, a corrected state estimate and likelihood score are computed. If no measurement 
falls within the gate, the track is not updated at that sample time. Once all tracks have been 
processed with the measurements for a given sample time, a track management routine 
begins. This routine performs two functions: track dropping and track association. The 
track dropping function checks likelihood scores and estimate update times. If these exceed 
thresholds, the tracks are dropped and will not be considered in future processing. The track 
association function compares the scores of all estimates that are candidates for the orginal 


two targets, and selects the estimates that correspond with the lowest scores. 


After looping through all measurement times, the simulation generates several 
output plots. True contact motion, noisy measurements, false measurements, predicted 
position estimates, and corrected position estimates are displayed in several figures. Also, a 
statistics matnx that is useful in tracker performance evaluation is displayed in the 
MATLAB command window. Items such as the number of tracks held, the number of 
clutter points generated, and the measurement and estimate error at each sample time can be 


reviewed. 


B. SIMULATION 1: FEATURE REJECTION TEST 
The first simulation demonstrates the ability of the feature rejection algorithm 


discussed in Section B of Chapter LI to discriminate between false measurements and target 


Zo 


measurements. Here, clutter frequency and PRF means are well separated from the two 
targets relative to the standard deviation of the sensor’s measurements. The parameters for 


the simulation are as follows: 
Simulation: 
number of measurements: 16 
clutter density: 5 false targets per 60nm X 60nm area 
Sensor: 
Position: 18.00N 113.00W 
O;; Snm 
Og. 0.5 degrees 
Ofreg. 9.001 GHz 
Oprf. 9.1 pps 
Mean Clutter Features: 
Setutter. 1.4 GHz 


PY cluner: 3.8 Pps 


Target 1: 





Initial Position: 22.00N 110.00W 


Heading: 145T 


24 


Speed: 13 Kts 
Tums: To 180T at 590 minutes in simulation 
Smean: 1.406 GHz 
PYfmean: 4.4 pps 
pelo 
Target 2: 
Initial Position: 22.00N 109.00W 
Heading: 120T 
Speed: 15 Kts 
Turns: none 
teea,: 1.594 GHz 
eijmean: 3-2 Pps 
q? x00 
Measurements: 14 cartesian measurements at about 90 minute intervals. Additionally, one 


Line of Bearing measurement will be obtained 30 minutes before the 5" cartesian 


measurement. A second line of bearing measurement will be obtained 45 minutes 


before the 10th cartesian measurement. 
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Filter Settings and Thresholds: 





Geographic gating: y = 6, corresponding to 95% probability for a chi-squared 
distributed random variable. 

Delete track score threshold: 48 

Delete track time-late threshold: 200 minutes 


Figure 2 shows the true track for targets one and two. Target one will follow the 


course change model, and target two will maintain a straight line track. 
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Figure 2. True Target Motion for Targets 1 and 2. 
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Figure 3 shows the true target position at sampling times, target and false 
measurements, and tracker position predictions (Zyy.;). The ellipses represent the 95% 
confidence regions around the predictions and are obtained from Pyy.; . These ellipses are 
related to the geographic gate used by the tracker. It is noted that many measurements are 
geographically feasible candidates, but it will be shown that, for this simulation, the feature 


rejection loop effectively removes most from consideration. 


Target Motion and Measurements at Sample Times 
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Figure 3. True Target Position and Measurements At Sample Times. 
The ellipses are obtained from the prediction covariance and are related to the 
geographic gate. 


Figure 4 shows true target position at sampling times, target measurements, and 


tracker corrected estimates (zyx). In this figure, the error ellipses represent the 95% 


Zi 


confidence region around the corrected estimate and are obtained from Py, . The longer 
Lines of Bearing correspond to target 2. The tracker is performing reasonably well with 


both targets. Estimate error increases during target 1’s turn but quickly diminishes. 


Target Motion, Measurements, and Estimates at Sample Times 
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Figure 4. Simulation 1 True Target Position, Target Measurements, and Corrected 
Estimates at Sample Times. 


Figure 5 contains the plot of Figure 4 for target 1 only. Likewise, Figure 6 1s the 


target 2 portion of figure 4. 
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Target 1 Motion, Measurements, and Estimates at Sample Times 
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Figure 5. Simulation 1 True Target Position, Target Measurements, and Corrected 


Estimates at Sample Times for Target 1. 
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Target 2 Motion, Measurements, and Estimates at Sample Times 
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Figure 6. Simulation 1 True Target Position, Target Measurements, and Corrected 


Estimates at Sample Times for Target 2. 
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Table 1 is a compilation of statistics from the simulation. The following is an 

explanation of the fields: 

index: measurement number 

report type: "1" indicates a cartesian measurement sensor report (Lat/Long). "0" 

indicates a line of bearing measurement sensor report. 

No of Clutter Msmts: The number of false measurements generated at that 
measurement time. 

No of Valid Tracks: The number of tracks being processed by the tracker. It 
includes the original two targets plus new tracks that were generated by the 
track splitting algorithm, less tracks that were deleted for exceeding score or 
maximum time without update thresholds. 

The following statistics are compiled for both targets: 

No w/ Valid Attributes: The number of measurements at the sampling time that 
have emitter frequency and PRF parameters close enough to the target's 
expected values to pass the feature rejection loop. 

No of Valid Measurements: The number of measurements at the sampling time to 
pass the feature rejection loop and fall within the range gate. If this number 
is "1", the track is updated with this measurement. If this number is greater 
than one, the measurement closest to the prediction will update the track, 
and remaining measurements will be used to generate new tracks (track 
split). If this number is "0", the track will not be updated at the sample time 


and output plots will contain the predicted, vice corrected, estimate. 


31 


Msmt error: Geometric distance in 60 nm (one degree) units between the target 
measurement and the true target position. 

Estimate error: Geometric distance in 60 nm (one degree) units between the 
corrected estimate and the true target position. 

Estimate swap: The track number of the candidate estimate with the lowest 
likelihood score. If this number is the same as the target number, the 
estimate for the target is the best available and no swapping will occur. If, 
however, this track number 1s different from the target number, the track 
numbers will be swapped so that the estimate for target one has the lowest 
likelihood score. 

No of Trk Candidates: The number of tracks that are candidates for this target. 
These include the original track plus all tracks that were generated by valid 
measurements (those passing the feature rejection loop and falling within the 
gate), less those tracks that were deleted for exceeding score or maximum 
time without update thresholds. 

Examination of the highlighted rows in Table 1 reveals the effectiveness of the 
feature rejection loop when target and clutter frequency and PRF means are well separated 
relative to the sensor’s standard deviations. Despite a large number of measurements not 
originating from the target, few have frequency and PRF measurements close enough to the 
target expected values to be considered valid. In the case of the simulation results presented 
here, typically only one measurement successfully passes this loop per target, presumably 


originating from the target being processed. 


index 

report type 

no of clutter msmts 

no of valid tracks 
target 1 

no of valid attributes 

no of valid msmts 

msmt error 

estimate error 

estimate swap 

no of trk candidates 
target 2 

no of valid attributes 

no of valid msmts 

msmit error 

estimate error 

estimate swap 

no of trk candidates 


index 

reporttype 

no of clutter msmts 

no of valid tracks 
target 1 

no of valid attributes 

no of valid msmis 

msmt error 

estimate error 

estimate swap 

no of trk candidates 
target 2 

no of valid attributes 


noofvalidmsmts 


msmt error 
estimate error 
estimate swap 

no of irk candidates 


1 2 3 4 5 6 7 
1 1 1 { 0 1 1 
0 0 ) 0 Ouumysonemed> 
2 2 2 2 2 2 2 
1 1 1 1 1 1 1 
J 1 1 oe 1 1 
0.1065 0.0537 0.1439 0.1056 NaN 0.0908 0.0176 
0.0900 0.0170 0.1175 0.1303 0.1725 0.0995 0.0302 
1 1 1 1 1 1 1 
1 1 1 1 1 1 1 
1 1 1 I 1 1 
{ 1 1 1 1 1 1 


0.0283 0.1311 0.0162 0.0448 NaN 0.0324 0.0693 


0.0529 0.0977 0.0399 0.0385 0.1040 0.0339 0.0441 
2 2 2 2 2 2 2 

1 1 1 1 1 1 1 

9 10 11 12 13 14 15 

GO pra O ane 0 eed 

2 2 2 2) a 2 2 

1 Te 3 1 1 1 1 

ies 1 mene | pam i 1 1 
0.0438 NaN 0.1330 0.0712 0.1047 0.1283 0.1503 
0.0846 0.1930 0.0808 0.0869 0.1089 0.1281 0.1492 
1 1 1 1 1 1 1 

1 1 1 1 1 1 1 

1 1 | ee 1 

fl — 1 1 1 1 
0.0539 NaN 0.1232 0.0918 0.1153 0.1139 0.1031 
0.0635 0.0758 0.1120 0.0588 0.0645 0.0704 0.1090 
2 2 2 2 D 2 2 

1 1 1 1 1 1 1 


Table 1. Statistics for Simulation 1: Feature Rejection Test 
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C. SIMULATION 2: GATING TEST 

Another case to consider is when the two targets have similar frequency and PRF 
means. This will demonstrate the effectiveness of the geographic gating in selecting the 
correct candidate measurement. This simulation is conducted with the same parameters as 


the previous simulation, with the exception of the target feature means: 


frean. 1.406 GHz 


PY} mean. 4.4 pps 


for both targets. The false measurement feature means remain unchanged, ensuring that the 
majority of these measurements will be rejected. 

Figure 7 shows true target position at sampling times, target measurements, and 
tracker corrected estimates (zx). As with figure 4, the error ellipses represent the 95% 
confidence region around the corrected estimate and are obtained from Py, . Plots of true 
target track and clutter measurements have been omitted since they are similar to Figures 2 
and 3. As expected, the tracker is performing reasonably well with both targets and 
producing results that are similar to the case when features are well separated. 

Analysis of the simulation results contained in Table 2 shows that, generally, the 
gating is effectively selecting the correct measurement. At each sample time, there are two 
mesurenionte that have passed through the feature rejection loop, presumably 
corresponding to targets 1 and 2. Gating usually selects the correct measurement for use in 


updating the estimate, preventing a track split. Occasionally however, more than one 
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Target Motion, Measurements, and Estimates at sample Times 
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Figure 7. Simulation 2 True Target Position, Target Measurements, and Corrected 
Estimates at Sample Times. 
measurement is in the geographic gate and the track is split. These instances are limited to 
the Line of Bearing measurements (indices 5 and 10), and should be expected for the 
geometry of the simulation. Referring to Figure 3, the line of bearing measurements to both 
targets fall within the ellipses for both predicted estimates, indicating a track split situation. 
The track is split by updating the prediction with both line of bearing measurements, one for 
each track. The generated tracks have state estimates that are very close to the updated 


track, making them nearly indistinguishable in figure 7. 
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index 

report type 

no of clutter msmts 

no of valid tracks 
target 1 

no of valid attributes 

no of valid msmts 

msmt error 

estimate error 

estimate swap 

no of trk candidates 
target 2 

no of valid attributes 

no of valid msmts 

msmt error 

estimate error 

estimate swap 

no of trk candidates 


index 

report type 

no of clutter msmts 

no of valid tracks 
target 1 


no of valid attributes | 


no of valid msmts 

msmt error 

estimate error 

estimate swap 

no of trk candidates 
target 2 


no of valid attributes 


no of valid msmts _ 
msmt error 
estimate error 
estimate swap 

no of trk candidates 


1 2 3 4 5 6 7/ 8 
1 1 1 1 0 1 1 1 

0 0 0 0 0 74 eG At 

p, D 2 2 4 4 4 4 

2 2 a5 2 2 2 oa 2 

1 1 1 1 2 1 14 1 
0.1584 0.0286 0.0336 0.1654 NaN 0.1298 0.0353 0.1480 
0.1174 0.1008 0.0581 0.1168 0.2270 0.1458 0.0280 0.1222 
1 1 1 1 1 1 1 1 

1 1 1 1 2 2 2 2 

>) ~~ Pew? PO ll. 2 ~p 2 
0.0638 0.1803 0.0583 0.1018 NaN 0.0745 0.0354 0.0529 
0.0341 0.1815 0.0528 0.0908 0.1547 0.0751 0.0189 0.0485 
2 2 2 2 2 2 2 2 

1 1 1 1 2 2 2 2 

9 10 11 12 13 14 15 16 

1 0 1 1 1 1 1 1 

42 42 74 115 45 43 40 A7 

4 8 8 8 8 8 8 8 

” oe 2 ao We 2 
—IF Fe Al ol! Jl eee 1 
0.0936 NaN 0.1133 0.1279 0.0443 0.0564 0.1033 0.0692 
0.1283 0.3668 0.4907 0.1643 0.0224 0.0620 0.0620 0.0634 
1 1 1 1 1 1 5 1 

2 4 4 4 4 4 4 4 
ea) ma 2 2 2 2 2 
0.1269 NaN 0.1548 0.1568 0.0887 0.0255 0.0426 0.2025 
0.0793 0.1603 0.1153 0.1263 0.0411 0.0387 0.0301 0.1650 
2 2 Z Dan ae 2 2 2 

z, 4 4 4 4 4 4 4 


Table 2. Statistics for Simulation 2: Gating Test 
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D. SIMULATION 3: SCORING TEST 
The tracker’s ability to select the best track candidate is examined in a simulation 
where the clutter’s feature means are relatively close (within 3 standard deviations) of a 


target’s means: 
Mean Clutter Features: 
ereres. 1.4 GHz 


P Wciutter: 3.8 Pps 


Mean Target 1 Features: 





Wmeane. 1.403 GHz 


PYfmean. 4.1 pps 


The tracker will have several false measurements successfully pass the feature rejection 
loop for target 1. Of these measurements, it 1s expected that some will fall within the 
geographic gate causing the track to be split. With successive measurements, track scoring 
is used to select the estimate that most closely conforms to the motion model and eliminate 
unlikeiy tracks. The number of measurements for this simulation is increased from 16 to 20 
to better show the track management capabilities. Mean features for target 2 are as 


specified in Simulation 1 and are expected to produce similar results. 


Figure 8 shows true target position at sampling times, target measurements, and 
tracker corrected estimates (zy). All tracks are shown here, including those generated as a 


result of splitting the original tracks. However, error ellipses are plotted only around the 
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estimates most likely to be targets 1 and 2 (1.e., those two candidate tracks with the lowest 
track scores). It is apparent that several new tracks have been generated for target 1, and 
that the track with the lowest score conforms to the true target track well. Figure 9 shows 
true target position at sampling times, target measurements, and tracker corrected estimates 


(Zxx) for the most likely track estimate only, illustrating this further. 


Target Motion, Measurements, and Estimates at Sample Times 
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Figure 8. Simulation 3 True Target Position, Target Measurements, and Corrected 
Estimates at Sample Times. 
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Target 1 Motion, Measurements, and Estimates at Sample Times 
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F igure 9. Simulation 3 True Target Position, Target Measurements, and Corrected 
Estimates at Sample Times. Best Estimate Only For Target 1. 
An understanding of the track management capability for this tracker can be obtained by 
examining Table 3. Multiple measurements possess valid features and pass through the 
feature rejection loop at many measurement times. Also, track splitting occurs when more 
than one measurement falls within the geographic gate: samples 9 and 12. These new 
tracks add to the list of track candidates for target 1. The number of track candidates 
increases from 1 to 3 at sample 9 (the first split), then peaks at 7 candidates at sample 13. It 
recedes to 2 candidates by the end of the simulation as a result of track dropping. By 


comparing the track score for all track candidates, the tracker identifies estimates that are 
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more likely to be target 1 than the currently held estimate on one occasion. At index 11, 
track 3 is selected and becomes the new track 1. The track shown in Figure 7 is the 


sequence of estimates with the lowest track score at each sample time. 


index 1 2 3 4 5 6 
report type 1 1 1 1 0 1 1 
no of clutter msmts 0 0 0 0 0 86 39 
no of valid tracks 2 2 2 2 2 2 Zz 
target 1 
no of valid attributes 1 1 0 1 1 P= m3 
no of valid msmts wail we ale eae | 1 al 1 
msmt error 0.0698 0.1495 0.13811 0.0244 NaN 0.0640 0.0719 
estimate error 0.03866 0.1360 0.2751 0.0138 0.0502 0.0518 0.0440 
no of trk candidates 1 - 1 1G 1 1 
target 2 
no of valid attributes 1 1 1 1 1 1 1 
no of valid msmts 1 1 1 1 1 1 1 
msmt error 0.0448 0.0685 0.0574 0.0721 NaN 0.1022 0.1017 
estimate error 0.0183 0.0648 0.0446 0.0500 0.0938 0.0941 0.0934 
estimate swap 2 2 2 Z Z 2 Z 
no of trk candidates 1 1 1 1 1 1 1 


Table 3. Statistics for Simulation 3: Scoring Test 
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index 

report type 

no of clutter msmts 

no of valid tracks 
target 1 

no of valid attributes 

no of valid msmts 

msmt error 

estimate error 

estimate swap 


no of trk candidates 


target 2 
no of valid attributes 
no of valid msmts 
msmt error 
estimate error 
estimate swap 
no of trk candidates 


index 

report type 

no of clutter msmts 

no of valid tracks 
target 1 

no of valid attributes 


no of valid msmts _ 


msmt error 

estimate error 

estimate swap 

no of trk candidates 
target 2 

no of valid attributes 

no of valid msmts 

msmt error 

estimate error 

estimate swap 

no of trk candidates 


{ 
0.0684 
0.0729 
eZ 
4 


S 10 

1 0 

43 43 

4 A 

o 5 

3 1 
0.0973 NaN 

0.3519 0.6189 

a 1 

3 3 

1 1 

1 1 
0.1509 NaN 

0.1233 0.1747 

2 2 

1 1 

16 17 

1 1 

46 46 

7 4 

8 2 

1 1 

0.0871 0.0535 

0.0205 0.0287 

1 1 

6 3 

1 1 

1 1 

0.0740 0.0996 

0.0490 0.0667 

Z Z 


4 


+ 


11 


0.1141 
0.1039 
$ 
4 


’ 
0.1714 
0.1177 
Bi 
1 


0.0624 
0:0572 
{ 
2 


' 
’ 
0.0962 
0.1037 
2 
1 


2 


48 


+ 


13 


46 


1 
0.1346 
0.1273 

; 

5 


’ 
’ 
0.0727 
0.0818 
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Table 3 (Continued). Statistics for Simulation 3: Scoring Test 
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’ 
, 
0.1237 
0.1055 
2 
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Repeated simulations with these initial parameters produce results that vary in detail 
but are generally the same in overall performance. The tracker consistently obtains a 
sequence of estimates that conforms with the true track reasonably well. Differences 
between simulations are found in the number and occurrence of track splitting. Also, some 
simulations have more estimate swaps than presented here, while others have less. This is 
expected due to the randomness of the clutter and target measurements. The track scoring 
and comparisons performed by the tracker have demonstrated considerable strength in 
separating sequences that arise from uniformly distributed clutter measurements from those 


that arise from Gaussian distributed target measurements. 


E. SIMULATION 4: COMBINED TEST 

Finally, the case is considered where false target features are broadly distributed 
over a range compared to the target features. This is the most realistic of the scenarios 
presented, since it models target tracking in an environment where many other contacts are 
operating with various emitters. A combination of the tracker’s capabilities will be required 
to maintain a reasonable track: feature rejection, gating, and track scoring. Furthermore, it 
will be assumed that the tracker does not know the feature mean for the false measurements 
and the alternate feature rejection technique discussed in Section B of Chapter LI will be 
used. For the purpose of generating false target feature measurements, a normal distribution 


is used with the following statistics: 


Mean Clutter Features: 
Setutter: 1.4 GHz 


PY cluuer 3-0 PPS 


Clutter Standard Deviation: 


Gireg. 1.000 GHz 


Oprf. 3.0 pps 


Target features are normally distributed as specified in Simulation 1: 
Target 1: 
tren |.406 GHz 
PYf mean: 4.4 pps 
Target 2: 
Reon le g94 GHz 


Die. PPS 


Target Standard Deviation: 





Ghreqg, 0.001 GHz 


Opry: 0.1 pps 
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Figure 10 shows true target position at sampling times, target measurements, and 
tracker corrected estimates (z,x). As before, all tracks are shown here with error ellipses 
plotted only around the estimates most likely to be targets 1 and 2. The tracker 1s 
performing reasonably well with both targets. Some new tracks have been generated 


around target 1, but the most likely state estimate conforms well with true target position. 


Target Motion, Measurements, and Estimates at Sample Times 
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Figure 10. Simulation 4 True Target Position, Target Measurements, and Corrected 
Estimates at Sample Times. 


Table 4 shows the output statistics for this simulation. It is noted that a few false 
measurements pass through the feature rejection loop frequently, but track splitting occurs 
rarely. There is only one track split apparent with target 1, and that occurs at sample 7. 
Target 2 does not exhibit any track splitting. The number of valid tracks increases from 2 to 
3 at sample 7, then from 3 to 4 at sample 11. Since there is only one valid measurement in 
the gate for both targets at this time, it suggests that the track generated at sample 7 was 
split. The number of valid tracks recedes to 3 at index 15, where it remains until the end of 
the simulation. It is also noted that none of the generated tracks ever achieve a lower score 


than the original two tracks and thus no estimate swapping occurs. 


index 1 2 3 4 S 6 ve 
report type 1 1 1 1 0 1 1 
no of clutter msmts @) 0 0 0 0 87 44 
no of valid tracks 2 2 2 2 2 2 S 
target 1 
no of valid attributes — «8 1 1 1 2 2 
no of valid msmts 1 1 1 1 1 1 2 
msmt error 0.0355 0.1014 0.1158 0.0805 NaN 0.0263 0.0098 
estimate error 0.0323 0.0898 0.1222 0.0978 0.1060 0.0288 0.0019 
estimate swap el 1 1 1 
no of trk candidates 1 1 1 = 1 1 eZ 
target 2 
no of valid attributes 1 1 1 1 1 3 2 
no of valid msmts 1 1 1 1 1 1 1 
msmt error 0.0421 0.0833 0.0262 0.0819 NaN 0.0804 0.0829 
estimate error 0.0374 0.0695 0.0372 0.0663 0.1019 0.0723 0.0725 
estimate swap 2 Z 2 2 2 z 2 
no of trk candidates 1 1 1 1 1 1 1 


Table 4. Statistics for Simulation 4: Combined Test 
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index 

report type 

no of clutter msmts 

no of valid tracks 
target 1 

no of valid attributes 

no of valid msmts 

msmt error 

estimate error 

estimate swap 


no of trk candidates 


target 2 
no of valid attributes 
no of valid msmts 
msmt error 
estimate error 
estimate swap 


no of trk candidates _ 


index 

report type 

no of clutter msmts 

no of valid tracks 
target 1 

no of valid attributes 

no of valid msmts 

msmt error 

estimate error 

estimate swap 

no of trk candidates 
target 2 


no of valid attributes - 


no of valid msmts 
msmt error 
estimate error 
estimate swap 

no of trk candidates 


1 
0.0436 
0.0329 
Z 
1 


9 10 11 
1 0 1 
42 42 81 
3 3 4 
2 “ay 1 
ie 1 1 
0.0953 NaN 0.0496 
0.0945 0.2575 0.0888 
1 1 1 
2 2 «83 
_— a 
1 1 1 
0.0371 NaN 0.0363 
0.0623 0.0786 0.0276 
2 2 2 
——e 
16 WW 18 
1 1 1 
44 42 50 
3 3 3 
1 ie 2 
i! a | 
0.0756 0.0542 0.0384 
0.0391 0.0600 0.0176 
1 1 i 
2 LB 
sr 5 
1 1 1 
0.0706 0.1265 0.0961 
0.0600 0.0749 0.0536 
2 2 2 
——1 1 


0.0846 
0.0564 
1 
3 


1 
—— 
0.1444 
0.1021 


0.1297 
0.0966 


0.1294 


0.0859 
2 
1 


13 


1 
0.0490 
0.0596 
1 
3 


3 
, 
0.1112 
0.1122 
2 


0.1203 
Oizge 
2 
1 


Table 4 (Continued). Statistics for Simulation 4: Combined Test 
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Just as with the previous simulations, repeated runs with these initial parameters 
produce results that vary in detail but are consistent in overall performance. Several runs 
were observed with no occurrence of track splitting. On the other hand, many simulations 
produced satisfactory results, but with track splitting and estimate swapping. The tracker’s 
proven capabilities in splitting, scoring, and swapping result in target estimates that conform 


reasonably well with truth. 


i SIMULATIONS WITH EXTERNAL MEASUREMENT DATA 

Much of the code contained in the Appendix is devoted to generating data that the 
tracker would receive from the sensor report. Further testing of the algorithm could be 
accomplished by stripping off this code and reading in collected sensor reports in the form 
of an external data file. With minor modifications, the tracker could be tested with position 
reports on two real tracks and several false measurements. Sensor reports must contain the 


following information: 


Time of Measurement 

Measurement Type: a flag for indicating Cartesian (Latitude and Longitude) or 
Line of Bearing processing. 

Position Measurement 

Error Ellipse: The sensor covariance can be derived from this as discussed in 


Chapter I. 
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Emitter Frequency 
Frequency Measurement Covariance 
Emitter PRF 


Emitter PRF Measurement Covariance 


Because the computer simulation in the Appendix is designed to evaluate the track 
maintenance capability of a track splitting approach, care must be taken in initializing the 
tracker. In the Appendix, the state estimates for Targets 1 and 2 are initialized by passing 
two noisy measurements of target position to the supporting m-file kal2init. The initial 
State estimate and covariance are then obtained as described in Chapter II. A similar 
approach could be taken when using external data by reading in two measurements for each 
target and passing them to kal2init. Another approach would be to simply provide the 
initial estimate as a priori knowledge and proceed directly to the main loop in the code. 
Similarly, emitter frequency and PRF means are provided as a priori knowledge in the 
Appendix code. When testing with external data, these a priori estimates must be 


consistent with the data for satisfactory performance. 
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V. CONCLUSION 


A. SUMMARY 


This research has developed a simulation to explore the use of a track splitting 
tracker in the maritime surveillance problem. The tracker uses Kalman Filtering to estimate 
target state from received measurements, which may be latitude and longitude, or line of 
bearing reports. Prior to filter update, measurements are filtered on the basis of features and 
geographic gates to eliminate unlikely measurement candidates. When more than one 
measurement falls within a track’s gate, the track is split to allow subsequent measurements 
to help determine validity. Track scoring with a likelihood function serves as the basis for 
determining which track estimate orginated from the modeled target's sequence of 
measurements. To prevent saturation of the tracker with false tracks, a track dropping 
routine eliminates the least likely estimates by considering track scores and time since last 


update. 


The effectiveness of measurement rejection on the basis of features was 
demonstrated in simulation by assigning feature means to the false measurements that were 
well separated from the true target means. Examination of the simulation results showed 
that few or no clutter measurements passed through this rejection loop, while the 
measurements originating from the target generally did. The Kalman Filter track estimates 
conformed well with the true target tracks. To investigate the effectiveness of the 


geographic gating, both targets were assigned the same feature means. Since the targets 
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were well separated geographically but had similar features, both target measurements 
successfully passed the feature rejection test and were then correctly paired with a track on 
the basis of the geographic gate. Finally, the ability of the tracker to select the best track 
estimate from the set of candidate tracks was tested by assigning similar feature means to 
the false measurements and one of the targets. During the simulation, the track was split 
several times as a result of false measurements passing both the feature rejection and 
geographic gate tests. The tracker produced an estimate of target track that conformed well 
with truth by selecting the lowest scored estimate at each measurement time. These 


simulations demonstrated several desirable features for a maritime surveillance tracker. 


B. FURTHER RESEARCH 


With minor modifications, the code in the Appendix can be tested with a real data 
set. Initially, testing should resemble the simulations presented here with two real targets 
and several false measurements at each measurement time. Target feature means must be 
known before hand, and target tracks must be initialized before clutter measurements are 
introduced. Such a test would confirm the viability of the track splitting tracker in a simple, 


real-world application. 


The algorithm can be significantly enhanced by adding a more sophisticated track 
initialization routine. Presently, the simulation tests the ability of the track splitting 
approach to maintain target tracks once obtained. Measurements that do not pass the 


feature test and fall within the geographic gate are discarded. An initialization routine 
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might establish tentative tracks with these measurements, which could then become regular 


tracks when scoring or number of update thresholds are met. 


A more accurate geographic model can be developed for this tracker. Here, 
distances are converted to geographic latitude and longitude with a flat earth model. Since 
sensors may operate over great distances in maritime surveillance, a spherical earth model 


should be used. 


St 


2 





APPENDIX 


This Appendix contains the track-splitting tracker simulation and supporting code. 


It was implemented in MATLAB® 5.1. 
Main Simulation Code: trackSplit.m 


YotrackSplit.m 

% This file simulates a multi-target, multi-sensor, track splitting tracker. Two 

% dimensional motion is simulated using a constant velocity discrete time target motion 
Yomodel for two contacts (Target 1 and Target 2). An instantaneous turn may be 
Zospecified by selecting a time to turn and entering that value for ’tumnTime” below. 
%Targets also have emitter frequency and prf attributes. | 

% Ameasurement of target position and attnbutes are obtained “about every 90 min”. 
% Both targets are measured at the same instant. Measurement noise is modeled as a zero 
Yomean Gaussian Random Variable with variance as set in the SENSOR part of the code. 
%The position measurements may be taken as Range and Bearing or Bearing only. 
%oMeasurements are taken as Range and Bearing except on measurement indices that 
JYohave been specified for Bearing only measurements. Also, an IF statement is in place 
Yofor use in “skipping” measurements (missed measurements). 

% Clutter measurements are also generated with target range and bearing 
Jomeasurements. Clutter measurements are modeled as uniformly distributed within a 
Yoregion around true track position. These measurements are assumed to have frequency 
Yoand prf attributes that are also modeled as Gaussian random variables. These statistics 
Yoare specified in the SENSOR part of the code. The clutter region is related to the 
Jeigenvalues of the estimate covariance matrix, P(k|k-1), which may extremely large for 
Zotracks in their infancy. For this reason, an IF statement controls when clutter 
Yomeasurements begin. 

% Clutter rejection is accomplished by comparing attributes of all measurements with 
Jothe track’s expected attributes. A weighted distance between clutter and target attribute 
Yomeans is calculated for each measurement and the ratio is formed. Measurements are 
Yorejected as clutter when the ratio exceeds a threshold, ’attThresh”. 

% From the remaining measurements, Geometric Association is accomplished by a 
Yocomparison is made by comparing expected target position (z(k|k-1) and measurement 
Joposition. The measurement that produces a track with the lowest log-likelihood score 
Jois selected for filter update. New Tracks are generated for all remaining measurements 
Jin the validation region (Chi-Square Distribution 95% Confidence Ellipse). A kalman 
%filter performs a correction step using the selected measurement, resulting in z(k|k). An 
Yoerror ellipse is obtained from the covariance of error matrix for plotting. 
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% ‘Track Dropping logic is included to control the number of tracks in the system. 

% After all tracks have been processed at a given measurement time, track scores and 
Yolast update time are compared to threshold values (maxScore and maxTime). Tracks 
Zothat exceed the threshold are dropped. 

% Output plots are: 


% True target motion 

% True target motion and measurements 

% True target motion, measurements, and corrected estimates (x(k|k)) with 
Zo error ellipse 

%o True target motion, measurements, and corrected estimates for target 1 
Jo True target motion, measurements, and corrected estimates for target 2 
Jo This simulation requires the following additional m-files: 

% kal2init.m {initializes the filter from 2 measurements } 

Jo elipa.m { obtains error ellipse from P matrix } 

% pole2cart.m {converts range and bearing measurements to 
% cartesian coords} 

%o get2d.m { generates cartesian position reports 

Jo (measurements) } 

% getLob.m { generates LOB position reports } 

% getClut.m generates “clutter” cartesian position reports } 
Jo getClutFeature.m { generates “‘clutter’’ features reports } 

clear 


%oToTo [0 To Mo To To To Yo Mo Mo %o % % Simulation and Sensor Initializationww~~AHHHD%A 
%Define number of measurements to be taken: 


nmeas = 16; Yonumber of measurements 

clutDens = 5; %odefines clutter density: 5 tgt/3600 sqr nm 
ZoSeNnsor 

slp = [-113.00;18.00]; Yosensor location in lon;lat deg 

sigmalr = 5/60; %2.5 nm of uncertainty in range 
sigmalb = .5*pi/180; %1/2 degrees uncertainty in brg 


slv = diag([sigmalr‘2;sigmalb*2]); Yocovaniance matnix for uncorrelated 
Yrange and brg msmts 


sigmaF = le-3; % 100 Hz uncertainty in f 

sigmaPrf = 0.1; %0.1 pps uncertainty in prf 
fClut = 1.400; %clutter feature statistical means 
priClut = 3.8: %Ghz and pps 


34 


sigmaFClut = 1; %clutter feature statistical STD 
sigmaPrfClut = 3; 
sigC = diag([1/sigmaFClut; 1/sigmaPrfClut]); %parameter “packaging” for use in 


sig = diag({1/sigmaF; 1/sigmaPrf]); Yoclutter rejection loop 

attClut = [fClut;prfClut]; 

c = sqrt(6); %oerror ellipsoids for 95% confidence 

gamma = 6; Joparameter for 95% threshold of chi squared dist 
attThresh = 10; Yothreshold for attribute discrimination 

maxScore = 48; % Track dropping parameters 


maxTime = 270; 
ToT To To To To Fo To To To To Vo Yo Vo Mo Vo ar get Initialization%%% %%%o%o% To ToT Vo To V0 % 
Yoinitialize the Target 1 state vector of positions and velocities 


hdg = 145*pi/180; % Target heading 

speed = 13/3600; %Target speed, kts 

x1 = [-110.00; Z%oInitial x posit (deg lon) 
speed*sin(hdg); ZYolnitial x speed (y/s) 
22.00; Yolnitial y posit (deg lat) 
speed*cos(hdg)]; JoInitial y speed (y/s) 

turnTime = 590; ZYoturn 590 minutes into simulation 


% Attribute Data for Target 1 


fTrue = 1.406; Yotrue emitter Freq (GHz) 

prfTrue = 4.4; Yotrue emitter prf (pps) 

attTrack(:,1) = [fTrue;prfTrue]; %’’ packaging” 

ZoInitialize the Target 2 state vector of positions and velocities 

hdg2 = 120*pi/180; %o Target heading 

speed2 = 15/3600; %oTarget speed, kts 

x12 = [-109.00; Zolnitial x posit (deg lon) 
speed2*sin(hdg2); Z%olnitial x speed (y/s) 
22.00; YoInitial y posit (deg lat) 
speed2*cos(hdg2)]; Z%olnitial y speed (y/s) 

turnTime2 = 1e6; Yono turn for tgt 2 


% Attribute Data for Target 2 


fTrue2 = 1.394; Yotrue emitter Freq (GHz) 
Pielme = or. Yotrue emitter prf (pps) 
attTrack(:,2) = [fTrue2;prfTrue2]; %” packaging” 


Jo To V0 Yoo ToT Yo Vo Yo Yo 70% MoTracker Parameter Initialization %%%%%%%M7%%V%% 
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% Define the H (measurement) matrix for cartesian measurement 
H = file@ 0 Q; 
0 0 ] 0]; 


%Define Target 1 covariance of plant noise matrix, Q 


qsaqr = le-8; %adjustable to account for 
Z%oplant model deviations 
del = 90: Yomean sample period 
Q = qsqr*[del*3/3 del*a2 0 0; 
del*2/2 del 0 O; 
0 0 del*3/3 delta, 
0 0 delp 22 del]; 


%Define Target 2 covariance of plant noise matrix, Q 


qsqr2 =le-8; 

del=o0: 

Q2 = qsqr2*[del43/3 delA2/2 0 OP 
dela2/2 de] 0 O; 
0 0 de]*3/3 del" 2/2 
0 0 dela22 del]; 


% Build the measurement time vector (msmts taken “‘about’” every 90 min) 
trbMeas = 90*ones(nmeas-3,1) + 4*randn(nmeas-3,1); Yorange/brg measmt 
JYointervals in minutes 


lob1l = rbMeas(5)-30; ZYoget a line of brg msmt 30 min 
%before 4" r/b msmt 
lob2 = rbMeas(9)-45; ZYoget a line of brg msmt 80 min 


%before 9" r/b msmt 


hit = [0; roMeas([1:4]); lob 1;rbMeas([5:8]);lob2;rbMeas([9:(nmeas-3)})]; 
%"*hit” is the vector of sample 
% intervals 
hit = round(hit); Yomeasmt intervals in integer mins 


To To To [0 [0% To Lo 0% MoM Track Initializationwr%wh%bbb%%%V%%o% VV %o% 
% initialize filter 1 with 2 measurements 
deltalnit = (90 + 4*randn(1)): Ysample interval for the 


56 


%two msmts used in init 


[xhat(:,1) phat xi] = kal2init(xi, H, Q, sigmalr, sigmalb, deltalInit, s1p); 
phatPack(:,1) = reshape(phat, 16,1); %”’phatPack”’ packages the 
Yocovariance matnx 
ZY into a column vector 


%initialize filter 2 with 2 measurements 
[xhat(:,2) phat2 xi2} = kal2init(xi2, H, Q2, sigmalr, sigmalb, deltaInit, s1p); 
phatPack(:,2) = reshape(phat2,16,1); 


YoInitialization for tracker and Output Matrices for track 1 


Zot = |]; %Cartesian msmts 

posout = []; Yotrue cartesian positions 

emer — |}: Yocartesian measurement error (from truth) 

time = 0; Yosample times 

zest = NaN*ones(4, 1); Yoposition estimate (H*xhat(k|k)) 

pa=ex; ZYostate initialization 

ind = 1 Yosample index 

oe = Ue %x-coord for LOB plotting 

yb]; %y-coord for LOB plotting 

xestEll = []; %x-coord for corrected estimate 
%o(zest) 95% ellipses 

yestEll = []; %oy-coord for corrected estimate 
%(zest) 95% ellipses 

xPredeil =|]: %x-coord for predicted estimate 
%(H*xPred) 95% ellipses 

yPredEll = []; Y%oy-coord for predicted estimate 
%(H*xPred) 95% ellipses 

esikimonr = |]; Yestimate error (from truth) 

vy = [I]; Y%ofeature measurement ([freq;prf]}) 

eclut@ut = []: Yoclutter cartesian positions 

yClutOut = []; %clutter freature measurements 

lam = zeros(2); Yotrack score and last update time 
%[score;updTime} 

nTracks=2; Yoinitial number of tracks 

valTrack = [1 2]; Y%oinitial vector of valid track numbers 

trkCand = [1 2]; J initial vector of track candidates 


%filter 2 initialization 


ZO0t2 =| 
posout2 = []; 
error2 = []; 


a7 


tine =): 

ZESt =| Ik 
Oy AW 
ee [I E 

ye2 = []; 

4 oe —al [| 

yb2 = []; 
Neste l=]: 
Vestal — |i 
xPredBil2 = fic 
VETeaiell2 = |]: 
esther = |}; 
y2=[]; 
clmt@nt =); 
yClutOut2 = []; 


ToT To To Yo To To To Yo Vo Yo %o Mo Yotracker simulation%% %% %o To %o To Go Yo To Yo Yo Yo Yo To To Yo Vo No 
% Loop through the simulation, generating target motion between 
Yosample times and taking measurements at hit times 


for 11 = 1:length(hit), 
il 
if (nTracks > 50) break,end; Yostop sim if capacity exceeds 50 tracks 


Gos RR eee Manel Motion / Filter Prediction ** *** F**K KE REA 
delta = hit(1); Yotime interval since last measurement: 


% Define the F Matrix (Transition Matrix) for discrete time 
% target motion with constant velocity 
P= eel delta QO 0; 

0 1 0 QO; 

0 0 ] delta; 

0 0 0 1]; 


Q = qsqr*[delta*3/3_ delta*2/2 0 0; 


delta*2/2 delta 0 Q; 
0 0 delta*3/3 delta*2/2; 
0 0 delta*2/2 delta]; 


%tilter prediction 

for tt(kNo = 1:nTracks 
phat = reshape(phatPack(:,trkNo),4,4); 
xpred(: ,trkKNo)=F*xhat(:,trkNo); 
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Ppred=F*phat*F’ + Q; 
PpredPack(:,trkNo) = reshape(Ppred, 16,1); 
end 


Y%obtain 3-sigma ellipses around estimate for display purposes 
if (11>5) 
Ppred=reshape(PpredPack(:,1),4,4); 
K = [Ppred(1,1) Ppred(1,3); 
Ppred(3,1) Ppred(3,3)]; 
[xell, yell, smaj, smin, theta] = elipa(K, c, xpred(1,1), xpred(3,1)); 
xPredEll = [xPredEIl xell’]; 
yPredEll = [yPredEIl yell’ ]; 


Ppred=reshape(PpredPack(:,2),4,4); 
K = [Ppred(1,1) Ppred(1,3); 
Ppred(3,1) Ppred(3,3)]; 
[xell, yell, smaj, smin, theta] = elipa(K, c, xpred(1,2), xpred(3,2)); 
xPredEIl2 = [xPredEI12 xell’}; 
yPredEII2 = [yPredEII2 yell’); 
end 


%target 1 motion step: x(k+1) = F*x(k) 
Yotarget changes course to 180 at sample time before tumTime 
if (time(length(time)) > turnTime) 
OD 10k 
x(4) = -speed; 
end 
X= Fx; 


Yotarget 2 motion step: x(k+1) = F*x(k) 
Y%otarget changes course to 180 at sample time before turnTime 
if (time(length(time)) > turnTime2) 
x2(2) = 0; 
x2(4) = -speed2; 
end 
y mas Se See 


Yotake measmts at hit times: 
%take measurement from current target state vector and append 
Yoto the true target position output matrix 


De 


time = [time time(length(time))+hit(ii)]; 


Ztrue = H*x; 
posout = [posout ztrue]; 
Zire? = ox 


posout2 =[posout2  ztrue2}; 
otek eee here eer nd Of Target Motion’ Batemerediemont cee 9. 4 


Op 7 26 28 A He EE EE 2 2K 2 2K KK EE EK NT egcurement Generation ** ***#** ** ### KK KK HK KK 


Yobtain noisy msmt from sensor 1 with associated covariance matrix 
Yin cartesian coords. 

%The “if” statement is used to simulate randomly missed measurements 
if (rand(1) < 1), 

sigR = sigmalr; 

sigB = sigmalb; 

sigma = [sigR; sigB]; Zopackaging 


%ofind difference in positions between sensor and tgt 
diff = ztrue - slp; 
diff2 = ztrue2 - slp; 


Yotarget | true bearing and range 
r = sqrt(diff(1)*2 + diff(2)*2); 
brg = atan2(diff(1),diff(2)); 


Yotarget 2 true bearing and range 
r2 = sqrt(diff2(1)*2 + diff2(2)2); 
brg2 = atan2(diff2(1),diff2(2)); 


Yotarget 1 attribute measurements 

fMeas = fTrue + randn(1)*sigmaF; 
prfMeas = prfTrue + randn(1)*sigmaPrf; 

y = [y [f{Meas;prfMeas]]; 

Yotarget 2 attribute measurements 

fMeas2 = fTrue2 + randn(1)*sigmaF; 
prfMeas2 = prfTrue2 + randn(1)*sigmaPrf; 


y2 = [y2 [fMeas2;prfMeas2]]; 
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Temi pede -((ii=—5)|(ii==—10)); Yospecify msmts to be LOB 
Jhere, 5" and 10” msmt 
Jowill be LOB, 
%all others cartesian 


if repTyp2d, %bearing and range msmt 
[z(:,1), R] = get2d(r, brg, sigma, slp); %tegt 1 
Yretums noisy msmt and 
Yomsmt covariance 
RPack(:,1) = reshape(R,4,1); 
ztilde = ztrue - z(:,1); 
disterror = sqrt(ztilde’*ztilde); %msmt error 


[z(:,2), R2] = get2d(r2, brg2, sigma, slp); ZYtgt 2 
Yoreturns noisy msmt and 
Ymsmt covariance 

RPack(:,2) = reshape(R2,4,1); 

Zulde2 = Zigue2? - ZZ): 

disterror2 = sqrt(ztilde2’*ztilde2);%msmt error 


%add the measurement to the measurement output matrix 
zout=[zout 2(:,1)]; 
ZOULZ = |(Zout2 Ae2)- 


%This block will produce clutter posit at attribute msmts 
Iia(iiees 
Ppred=reshape(PpredPack(:,1),4,4); 
S1 = H*Ppred*H’ +R; 
[clut, nClutter1] = getClut(S1, ztrue, clutDens, gamma); 


Yobtains clutter posit msmts 

Yuniformly dist around tgt 1 

yClut = getClutFeature(nClutterl, fClut, prfClut, sigmaFClut, 
sigmaPrfClut); %obtains clutter feature msmts 


clutOut = [clutOut clut]; 
yClutOut = [yClutOut yClut]; 


Ppred=reshape(PpredPack(:,2),4,4); 
S2 = H*Ppred*H’ + R2; 
[clut2, nClutter2] = getClut(S2, ztrue2, clutDens, gamma); 
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yClut2 = getClutFeature(nClutter2, fClut, prfClut, sigmaFClut, 
sigmaPrfClut); 


clutOut2 = [clutOut2 clut2]; 
yClutOut2 = [yClutOut2 yClut2]; 


else You<5 ==> do not generate clutter msmts 
nClutterl = 0; 
nClutter2 = Q; 
Clliicce tpl 
wis = ||] 
yClut = []; 
ye lne — |]: 
end Yoclutter msmts for 2d case 


else Zolob msmt 
[z(:,1), xLob, yLob] = getLob(r, brg, sigB, s1p); 
Yotgt 1 noisy bearing, 
Yoand coordinates for 


%oplotting LOB 


[z(:,2), xLob2, yLob2] = getLob(r2, brg2, sigB, slp); %tgt 2 


Yoadd the LOB measurement to the output LOB measurement matrix 


xb = [xb xLob’]; 

yb = [yb yLob’]; 

xb2 = [xb2 xLob2’]; yb2 = [yb2 yLob2’]; 

ztilde = nan; Yomsmt distance errors not computed 


% for LOB 
disterror = nan; 
Ztilde2 = nan; 
disterror2 = nan; 
theta = brg; 


end §%measurement generation 
O76 er eRe ee nok oh ee of Measurement Generation ®* ** *** #* #4 # # EK 


J%add the measurement distance error to measurement distance error matrix 
error= [error  disterror]; 
error2 = [error2 disterror2]; 
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fp 76 7 28 Ae Baia 7 2 2K 2 2 2K EET a Processing / Filter CGomection=~* 2°. 4. 
% Target Clutter Rejection: Select correct measurement from vector of 
Yomeasurements by examining attributes nClutter = nClutter1+nClutter2; 
pos Vec = [z(:,1); z(:,2); clut({1:2*nClutter1])’;clut2({1:2*nClutter2])’]; 

Yovector of msmt positions 
Yonote that the correct msmt 
%for target 1 is the 1* 
Ymsmt in posVec. Target 
%2 is second. All others 
Y%are clutter 


attVec = [fMeas: prfMeas; fMeas2; prf{Meas?2; reshape(yClut,2*nClutter 1,1); 


reshape(yClut2,2*nClutter2, 1)]; Yovector of msmt attributes 
for trkNo = valTrack Y%oloops through all tracks 
nVal =0; % initialize valid msmt counter 


Ppred = reshape(PpredPack(:,trkNo),4,4); 
%extract Ppred for track being 
Yoprocessed from packed matrix 
for | = 1:(nClutter+2), 
att = [att Vec(2*]-1);attVec(2*1)]; 
dT = (att - attTrack(:,trkNo))’ *sigT*(att - attTrack(:,trkNo)); 
dC = (att - attClut)’*sigC*(att - attClut); 
dTdC(1) = dT/dC; 
end %this loop compares the distance of 
%a msmts features to the tgt and 
%clutter means. dTdC is a ratio of 
%these distances. 


valAtt = find(dTdC < gamma); YvalAtt holds the index of all msmts 
%that fall below the threshold. 
%These msmts have attributes that 
Y%are valid for this tgt 


alhdGaiir %reset for next 
%Geometric Track Association and Filter Update: associate observation that 


% gives lowest score for this track. 
% Track scoring determined by log-likelihood function 
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YoCASE 1: 2 Dimensional reports using standard filter 
YoCASE 2: LOB report using EKF 


xhatSplit(:,1)=[lam(1,trkNo)+6; xpred(:,trkNo); attTrack(2*trkNo-1); 
attTrack(2*trkNo)]; 
YoxhatSplit is a holding vector that will contain 
Zothe score, state vector, and attributes that result 
Z%ofrom all measurement associations with the track 
Ybeing processed. Here, it is initialized with the 
Yotracks score, predicted state, and attributes. 


phatS plitPack(:,1) = [lam(1,trkNo)+6; PpredPack(:,trkNo)]; 
%likewise, phatSplitPack holds the covariances 
Z%oresulting from an association 


if repTyp2d %2d (cartesian) report 


R = reshape(RPack(:,trkNo),2,2); %extract sensor covariance 
S = H*Ppred*H’ + R; %S is the innovation covar 


for 1 = 1:length(valAtt) 
ZObs(:,1) = [pos Vec(2*valAtt(1)-1);pos Vec(2*valAtt(1))]; 
ZYoobtain a valid (passed the attribute test) msmt 


ztil = zObs(:,]) - H*xpred(:,trkNo); %the innovation 


lamSplit = ztil’*inv(S)*ztil; Yonorm of innovation 
Yosquared, a distance 
Yomeasure 

if (lamSplit < gamma) %if the dist is inside gate 


lam(2,trkNo) = time(length(time)); 
Yoindicate time of update in 
Joscore matrix 

nVal =nVal-+ 1; JY indicate # of valid associat 


G = Ppred*H’*inv(S); %filter gain 
xhatSplit(:,n Val) = (lam(1,trkKNo)+lamSplit; 
(xpred(:,trkKNo) + G*ztil); 
attVec(2*valAtt(1)-1); 
attVec(2*valAtt(1))]; 
Zonow, xhatSplit holds new score and 
Yocorrected est 
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phatSplit = (eye(size(Ppred))-G*H)* Ppred 
*(eye(size(Ppred))-G*H)’+G*R*G’ ; 
phatSplitPack(:,nVal) = [lam(1,trkKNo)+lamSplit; 
reshape(phatSplit,16,1)]; 
end %valid msmt assoc 


end %72d report association 


else 


%LOB geom association and update 
R = sigB; Yosensor covar 
x1 = xpred(1,trkNo)-s1p(1); 


x3 = xpred(3,trkNo)-s1p(2); 
Cen =< 2 + x572: 


h={[-x3/den 0 x1/den 0]; 

%\inearized msmt matrix 
zhat = atan2(x1,x3); Yomsmt est 
S =h*Ppred*h’ + R; Yoinnovation covar 


for 1 = 1:length(val Att) 
ZObs(:,1) = [pos Vec(2*valAtt(1)-1); pos Vec(2*valAtt(1))]; 
%obtain a valid observation 


ztil = zObs(1,]) - zhat; %innovation 
lamSplit = ztil’*inv(S)*ztil; © %norm of innovation 
ZJosquared 


if (lamSplit < 4) Zoif inside gate... 
lam(2,trkKNo) = time(length(time)); 
Yoindicate time of update 
%in score matrix 
nVal = nVal + 1;%indicate # of valid assoc 


G = Ppred*h’*inv(S); %gain 
xhatSplit(:,nVal) = [lam(1,trkKNo)+lamSplit; 
(xpred(:,trkKNo) + G*ztil); 
att Vec(2* valAtt(1)-1); 
attVec(2*valAtt(1))]; 
Zonow, xhatSplit holds new score 
Yand corrected est 
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phatSplit = (eye(size(Ppred))-G*h) * Ppred * 
(eye(size(Ppred))-G*h)’ + G*R*G@’; 


phatSplitPack(:,nVal) = (lam(1,trkKNo)+ 
lamSplit; 
reshape(phatSplit,16,1)]; 
end 
end 
end %LOB geometric assoc and update 


Jo%% 
if trkKNo==1 %this IF block obtains the no 
valAttl = length(valAtt); Yoot msmts that pass the 
% feature clutter rej- 
nVall =nVal; %ection and geographic gates 
%for tracks 1 and 2. 
elseif trkNo== 
%o These values will be passed 
%to the matnx STAT 
valAtt2 = length(valAtt); %for examination in the 
nVal2 = nVal; %command window after 
%simulation is concluded 
end 


Yoafter simulation, this info can be reviewed in the command 

Ywindow. For trk 1, the index “1” should be in the valAtt vector. If 
Yonot, att Thresh may be too high. Likewise, index “2” should be included 
%in valAtt for trkNo 2. 


xhatSplit = sortrows(xhatSplit’);%sorts validated estimates into 

attSplit = xhatSplit(:,[6:7])’; ascending order of likelihood 
%function. attSplit holds the 
Yattnibute msmt 


phatSplitPack = sortrows(phatSplitPack’); %sorting by likelihood 
%here keeps covariances 
Z%indexed with estimates 


xhat(:,trkNo) = xhatSplit(1,[2:5])’; 

phatPack(:,trkKNo) = phatSplitPack(1,[2:17])’; 

lam(1,trkKNo) = xhatSplit(1, 1); 
%\lowest score is chosen for tgt track. 
%Estimate, Covar, and score are 
Youpdated. 
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for 1 = 2:nVal Yoall remaining associations 
d = size(xhat); Yinitialize new tracks. 
xhat(:,d(2)+1)=xhatSplit(1,[2:5])’;  %...new trk state est 
phatPack(:,d(2)+1)=phatS plitPack(1,[2:17])’; %...new trk 


%est covar 
lam(1,d(2)+1)=xhatSplitd,1); %...new trk score 
lam(2,d(2)+1)=time(length(time));  %...new trk init time 
att Track(:,d(2)+1)=attSplit(:,1); %...new trk att 
RPack(:,d(2)+1) = RPack(:,trkKNo); | -%...new trk sensor 

Yocovar 
zest(2*d(2)+1,:) = nan; %...new track in estimate 
zest(2*d(2)+2,:) = nan; Youtput matrix 
trkCand(d(2)+1) = trkCand(trkNo); 
valTrack = [valTrack d(2)+1]; 

end 
xhatSphit = []; ZYoretums this to a column vector 
phatSplitPack = []; 

end Z%trk processing loop 


%%%%Track Management Algorithm: 2 functions%%% 
%1) Drop Tracks that exceed thresholds (score and last update). 
%2) Select candidate with lowest score to be target estimate 
% After all tracks have been processed for a sample time... 
Zotrack dropping... 
delTrk = find((lam(1,:)>maxScore) | (time(length(time))-lam(2,:))>maxTime); 
%form a vector of track numbers for track scores 
Z%that exceed maxScore or update times that exceed 
Jomax Time 
if delTrk 
for | = delTrk 
valTrack(find(valTrack==])) = []; 
%deletes track from vector of 
end Yovalid tracks 
end 


xhat(:,delTrk)=nan; 
phatPack(:,delTrk)=nan; Y%Not-a-number indicates “track deleted” 
lam(:,delTrk)=nan; 


Zobest track selection... 

Zothis gets the track number of the candidate w/ lowest score 
new 1 = floor(finddam==min(lam(1,find(trkCand==1))))/2)+1; 
new2 = floor(find(lam==min(lam(1,find(trkCand==2))))/2)+1; 
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Yothese lines swap the state vector and covariance matrices 
Zoto put the “new” (lower score) vectors into trkNo | position. 
oldxhat=xhat(:,1); 

oldphat=phatPack(:,1); 

oldlam=lam(:,1); 


xhat(:,1) = xhat(:,new1); 
phatPack(:,1) = phatPack(:,new1); 
lam(:,1) = lam(:,new1); 


xhat(:,new1) = oldxhat; 
phatPack(:,new1) = oldphat; 
lam(:,new 1) = oldlam; 


%... likewise for trkNo 2 
oldxhat=xhat(:,2); 
oldphat=phatPack(:,2); 
oldlam=lam(:,2); 


xhat(:,2) = xhat(:,new2); 
phatPack(:,2) = phatPack(:,new2); 
lam(:,2) = lam(:,new2); 


xhat(:,new2) = oldxhat; 
phatPack(:,new2) = oldphat; 
lam(:,new2) = oldlam; 


else Zomissed msmt 
error = [error nan]; 
end 


Oo 78 75 28 EA 2K EK A A KE ACK EEF yr output Matrices* * * *# # KKK KKK HK HK KK 


d = size(xhat); 
nTracks = d(2); %o# columns in xhat is # of tracks held 


Yoadd estimate to the matrix of estimates 
for l=1:nTracks 
Ze( 1) = A xmalicl): 


end 


zest = [zest reshape(ze,2*1,1)]; 
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%obtain target | error ellipse of estimate 
phat = reshape(phatPack(:,1),4,4); 
K = [phat(1,1) phat(1,3); 
phat(3,1) phat(3,3)]; 
[xell, yell, smaj, smin, theta] = elipa(K, c, xhat(1,1), xhat(3,1)); 
mest —||xestill xell’ |; 
yestEll = [yestEll yell’); 
ze=[zest(1,length(zest(1,:))); zest(2,length(zest(2,:)))] ; 
estError = [estError sqrt((ztrue-ze)’*(ztrue-ze))]; 


%obtain target 2 error ellipse of estimate 
phat2 = reshape(phatPack(:,2),4,4); 
K = [phat2(1,1) phat2(1,3); 
phat2(3,1) phat2(3,3)]; 
[xell, yell, smaj, smin, theta] = elipa(K, c, xhat(1,2), xhat(3,2)); 
xestEll2 = [xestEII2 xell’]; 
yestEll2 = [yestEll2 yell’); 
ze=[zest(3,length(zest(3,:))); zest(4,length(zest(4,:)))] ; 
estError2 = [estError2 sqrt((ztrue2-ze)’ *(ztrue2-ze))]; 


Z%compile statistics for output 


stat(1 ,11)=11; %the matrix STAT will be displayed in the 
stat(2,11)=rep Tl yp2d; Yocommand window for post-simulation 
stat(3,11)=nClutter; Yoperformance analysis 
stat(4,i1)=length(valTrack); 

stat(5,11)=nan; 


stat(6,11)=valAtt]; 

stat(7,11)=n Vall; 

stat(8,11)=disterror; 
stat(9,11)=estError(length(estError)); 
stat(10,11)=new1; 
stat(11,11)=length(find(trkCand(valTrack)==1)); 
stat(12,11)=nan; 

stat(13,11)=valAtt2; 

stat(14,11)=n Val2; 

stat(15,11)=disterror2; 
stat(16,11)=estError2(length(estError2)); 
stat(17,11)=new2; 
stat(18,11)=length(find(trkCand(valTrack)==2)); 


end ZYotime step 
ToT To To Mo To To Vo Yo Yo Vo To Fo VoVoO ut put Results%%%% %o Go Vo No To To To Vo To Yo Yo Yo 
time(1) = []; 
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%Plot the results 
nEst = nmeas; 
nmeas = max(size(zout)); 


% Target motion and true position at measurement times 

figure(1) 

axis equal 
plot(posout(1,:),posout(2,:),’-",posout2(1,:),posout2(2,:),’-‘,slp(1),s1p(2),’ X’);gnid; 
title(“True Target Motion’); 

xlabel(‘Deg Lon’);ylabel( “Deg Lat’) 


%Target motion with noisy measurements 

figure(2) 

nbrg = min(size(xb)); 

nPred = min(size(xPredEl])); 

plot(posout(1,:),posout(2,:),’o’,zout(1,:),zout(2,:),’*’ .xb(:,[1:nbrg]), yb(:,[1:nbrg]), 
xb2(:,[1:nbrg]),yb2(:,[1:nbrg]),clutOut(1,:),clutOut(2,:),’*’,s1p(1),slp(2),’X’); 

hold on 

plot(posout2(1,:),posout2(2,:),’o’ ,zout2(1,:),zout2(2,:),’*’ ,clutOut2(1,:),clutOut2(2,:),’*’, 
slp(1),s1p(2),’X’,xPredEIl(:,[1:nPred]),yPredEIl(:,[1:nPred]), 
xPredElI2(:,[1:nPred]),yPredEIl2(:,[1:nPred]),slp(1),sIp(2),’X’);grid; 

axis equal 

title(“Target Motion and Measurements at Sample Times’); 

xlabel(“Deg Lon’);ylabel(“Deg Lat’) 

hold off 


% Target motion with Filter Estimates 

figure(3) 

plot(posout(1,:),posout(2,:),’o’ ,zout(1,:),zout(2,:),’*’,xb(:,[1:nbrg]),yb(:,[1:nbrg]), 
xb2(:,[1:nbrg]),yb2(:,[1:nbrg]),zest(2*[1:nTracks]-1,:),zest(2*[1:nTracks],:),’x’, 
xestEll(:,[1:nEst]),yestEll(:,[1:nEst]),s1p(1),s 1p(2),’X’);grid; 

hold on 

plot(posout2(1,:),posout2(2,:),’o’ ,zout2(1,:),zout2(2,:),’*’, 
xestEll2(:,[1:nEst]),yestEll2(:,[1:nEst]),s1p(1),s1p(2),’X’ );grid; 

axis equal 

title(“Target Motion, Measurements, and Estimates at Sample Times’); 

xlabel(“Deg Lon’);ylabel(“Deg Lat’) 

hold off 


% Target 1 motion with Filter Estimates 

figure(4) 

plot(posout(1,:),posout(2,:),’o’ ,zout(1,:),zout(2,:),’*’ ,xb(:,[1:nbrg]), yb(:,[1:nbrg]), 
zest(1,:),zest(2,:),’x’ ,xestEll(:,[1:nEst]), yestEll(:,[1:nEst]),slp(1),s1p(2),’X’);gnd; 
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axis equal 
title “Target 1 Motion, Measurements, and Estimates at Sample Times’); 
xlabel(“Deg Lon’);ylabel(“‘Deg Lat’) 


%o Target 2 motion with Filter Estimates 

figure(5) 

plot(posout2(1,:),posout2(2,:),’o’ ,zout2(1,:),zout2(2,:),’*’ ,xb2¢,[1:nbrg]),yb2(:,[1:nbrg]),ze 
st(3,:),zest(4,:),’x’ xestEll2(:,[1:nEst]), yestEll2(:,[1:nEst]),s1p(1),s1 p(2),’X’); grid; 

axis equal 

title(‘Target 2 Motion, Measurements, and Estimates at Sample Times’ ); 

xlabel(“Deg Lon’);ylabel( “Deg Lat’) 


Stat 


kal2init.m 


function [xhat, phat, xiNew] = kal2init(xi, H, Q, sigR, sigB, delta, sensPos) 
%initializes the kalman filter based on 2 measurements 

% 

%The 2 "true" positions are xi and F*xi. 

%oThese are converted to noisy measurements using pol2cart. The associated 
JYocartesian covariance matrices are also obtained. 


%From these noisy measurements and covariances, the initial state and covariance 
%ofor the filter is estimated. 


%INPUT parameters: 

Jo XI: true target state at time zero 

Yo tele Measurement matrix 

%o Q: Plant Noise 

Yo sigR: Sensor Range Inaccuracy (std dev) 
Jo sigB: Sensor Bearing Inaccuracy (std dev) 
Jo delta: Sample Interval 

Yo sensPos: Sensor Position 

%Y OUTPUT parameters 

Yo xhat: Initial State Estimate 

% phat: Initial State Estimate Covariance 

Jo xiNew: new true target state (x1 at time delta) 


% Supporting m-files: 
% pole2cart.m 
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% Define the F Matrix (Transition Matrix) for discrete time 
% target motion with constant velocity 
F=[1 delta O 0; 

0 1 0 0; 

0 0 1 delta; 

0 0 0 iN 


(z1 R1] = pole2cart(H*xi, sensPos, sigR, sigB); 
xiNew = F*xi; 
[z2 R2] = pole2cart(H*xiNew, sensPos, sigR, sigB); 


72 = 77 71} 
A= 0 0 0; 
1/delta 0 -1/delta0; 
0 1 0 0; 
0 1/delta O -1/delta]}; 
mv’  eeaimv(): 
xhat =A*z: 


phat = A*[R2 zeros(2,2);zeros(2,2) R1+H*invf*Q*invf*H7]*A* 


pole2cart.m 
function [z, R, disterror] = pole2cart(ztrue, sp, sigmar, sigmab); 


Yobtains a noisy polar measurement to target, adds measurement noise, 
Yoand converts measurement and covariance matrix to cartesian 


Zoinputs: 

Jo ztrue = [xt; yt] true posit of target in cartesian 

Jo Sp = [xs; ys] posit of sensor in cartesian 

Jo sigmar = scalar range msmt standard deviation 

Jo sigmab = scalar bearing msmt standard deviation 
Yooutputs: 

Jo z= [xm; ym] noisy Measurement in cartesian 

Jo R=2x2 covariance matrix in cartesianj 

% disterror = sqrt((ztrue - z)*(ztrue - z)’) 


1? 


sigma = [sigmar; sigmab]; 


Yofind difference in positions between sensor and tgt 
diff = ztrue - sp; 


Yotrue bearing and range 
r= sqrt(diff(1)*2 + diff(2)2); 
brg = atan2(diff(1),diff(2)); 
Zplr = [r; brg]; 


%add noise to the measurement 
zmeas = zplr + randn(size(sigma)).*sigma; 


Yconvert the noisy measurement to cartesian 
zZ = [zmeas(1)*sin(zmeas(2)); 
zmeas(1)*cos(zmeas(2))]; 
Z—= Sp +Z; 


%convert covariance matrix to cartesian 
sv = diag(sigma.”*2); 
M = [sin(zmeas(2)) -zmeas(1)*cos(zmeas(2)); 
cos(zmeas(2)) | zmeas(1)*sin(zmeas(2))]; 
R = M*sv*M’; 
%calculate distance error 


Ztilde = ztrue - Z; 
disterror = sqrt(ztilde™ztilde); 


elipa.m 
function [xout, yout, smaj, smin, theta] = elipa(PK, c, xt, yt) 


Jocalculates error ellipsoids given error covariance 
Yoand estimate position 


%INPUTS: 

Jo PK is covariance matrix 

Jo c is conficence region (c=sqrt(6) for 95%) 
Jo Xt is ellipse center x coord 

Jo yt is ellipse center y coord 


ie. 


%OUTPUTS: 


Yo Xout 1s a vector containing x coord points of the ellipse 
Jo yout is a vector containing y coord points of the ellipse 
Jo smaj is the length of the semi-major axis 

Jo smin is the length of the semi-minor axis 

% theta is the orientation of the semi-major axis 


% to plot the ellipse, use: plot(xout, yout) 
Yadapted from Stephen L. Spehn’s Errellip.m (15 Nov 89) by Mark Olson 


%get eigenvalues (lam) and eigenvectors(V) 
[V,lam] = eig(PK); 

sigx = sqrt(lam(1,1)); 

sigy = sqrt(lam(2,2)); 


Yoparameterized ellipse 
t = 0:2*pi/100:2*pi; 

X = sigx*c*cos(t); 

y = sigy*c*sin(t); 


%translate to eigenvectors space and center at tgt posit 
xout = x* V(1,1) + y* V(1,2) + xt; 
yout = x*V(2,1) + y* V(2,2) + yt; 


Yoreport semimajor and semiminor axes lengths and onentation 

smaj = 2*c*max([sigx sigy]); 

smin = 2*c*min([sigx sigy]); 

theta = atan2(V(2,1),V(1,1)) + pi/2*(sigy > sigx); %sigy > sigxk ==> y 1S smaj 


get2d.m 
function [z, R] = get2d(r, brg, sigma, sensPos) 


%this function generates a noisy cartesian measurement from 
Yotrue bearing and range. Polar measurement noise is added to 
Yotrue polar position, then rotated into a cartesian coordinate 
Yoreference. The measurement covariance is also rotated into a 
Yocartesian coordinate reference. 


ZINPUTS: 
Zo I: true range to target 
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% brg: true bearing to target 


Yo sigma: column vector of sensor std dev [sigB;sigR] 
Jo sensPos: sensor position in cartesian coordinates 

% OUTPUTS: 

Yo Zs Cartesian measurement 

Yo R: cartesian measurement covariance 


Y%obtains cartesian measurement true range and brg msmt 
Zplr = [r; bre]; 


Yoadd noise to the measurement 
zmeas = zplr + randn(size(sigma)).*sigma; 


Yoconvert the noisy measurement to cartesian for plotting 
z = [zmeas(1)*sin(zmeas(2)); 


zmeas(1)*cos(zmeas(2))]; 
Z = sensPos + Z; 


%obtain target 1 covariance of error in cartesian 
M = [sin(zmeas(2)) -zmeas(1)*cos(zmeas(2)); 
cos(zmeas(2)) zmeas(1)*sin(zmeas(2))]; 


R = M*diag(sigma.42)*M’; 


getLob.m 
function [zbrg ,xLob, yLob] = getLob(r, brg, sigB, sp); 


Yothis function generates a noisy measurement bearing and output 
%vectors for plotting from true target position. 


%INPUTS 

% iP: true target range 

% brg: true target bearing 

%o sigB: sensor bearing measurement std dev 
% slp: sensor position in cartesian coord 

% OUTPUTS 

% zbrg: measured bearing 


a5 


Jo xLob: 
Jo yLob: 


Line of Bearing x-coords for plotting 
Line of Bearing y-coords for plotting 


Zoto plot the Line of Bearing, use 
% plot(xLob,yLob) 


%tind Line of Bearing Msmt to target 


%add noise to the measurement 
zbrg = brg + randn*sigB; 


Yoconvert the noisy measurement to cartesian for plotting 
xLob = linspace(3,(r+1)*sin(brg)); 
yLob = xLob/tan(zbrg); 
XxLob = slp(1)*ones(size(xLob)) + xLob; 
yLob = sl p(2)*ones(size(xLob)) + yLob; 


getClut.m 


function [clut, nClutter] = getClut(S, ztrue, clutDens, gamma); 


%this function generates uniformly distributed clutter position measurement 
Yin a region around the true target position. (Bar Shalom and Li, 1995) 


YINPUTS: 

% S: 

% ztrue: 

% clutDens: 
%o gamma: 
% 

YOUTPUTS: 

% clut: 

Yo 


% nClutter: 


Innovation Covariance S = H*P*H’+R 

true cartesian position of target, [x-coord;y-coord] 

Clutter Density, tgts/area (here, tgts/3600 sq nm) 
chi-squared dist threshold corresponding to validation region 
size. (here, gamma=saqrt(6) for 95% probability ellipse 


matrix of clutter positions clut = [x-coord of all clutter msmts; 
y-coord of all clutter msmts] 
number of clutter msmts generated 


Z%obtain clutter posit measmts 


Ygenerate target | clutter measurements: uniformly distributed in a 
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Yowindow around the true target position: 


Yodetermine the area of the validation region (ref: YBS MMT sec 3.4.6) 
Vk = pi*gamma*sqrt(det(S)); 


nClutter = round(10*Vk*clutDens+1);%Number of false msmsts in a square 
Yo around truth 


xClutter = sqrt(10* Vk)*(rand(1,nClutter) - .5*ones(1,nClutter)) + ztrue(1); 
yClutter = sqrt(10* Vk)*(rand(1 ,nClutter) - .5*ones(1,nClutter)) + ztrue(2); 


clut = [xClutter; yClutter]; 


getClutFeature.m 
function yClut = getClutFeature(nClutter, fClut, prfClut, sigmaFClut, sigmaPrfClut); 


Zthis file obtains the feature (emitter freq and prf) measurements for 
Zothe false targets. Normal dist 1s assumed. 


YINPUTS: 

% nClutter: the number of clutter targets 

% fClut: mean emitter frequency for clutter targets 

% prfClut: mean emitter prf for clutter targets 

% sigmaFClut: std dev of emitter freq for clutter targets 

% sigmaPrfClut: std dev of emitter prf for clutter targets 
%OUTPUT: 

% yClut: matrix of feature measurement column vectors 
% [emitter freq for all targets; 

%o emitter prf for all targets] 


fMeasClut = fClut*ones(1,nClutter) + randn(1,nClutter)*sigmaFClut; 
prf{MeasClut = prfClut*ones(1,nClutter) + randn(1,nClutter)*sigmaPrfClut; 


yClut = [fMeasClut;prfMeasClut]; 


ag. 
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